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Abstract 

This note details the development of a discrete-time diffusion process to approximate the 
midnight customer count process in a -Mper/Geo 2 timeScaie/-V system. We prove a limit theorem 
that supports this diffusion approximation, and discuss two methods to compute the stationary 
distribution of this discrete-time diffusion process. 


The Mper/Geo 2 timeScaie/-^ is siiigle-pool queueing system with a periodic Poisson arrival process 
and a two-time-scale service time feature. This queueing system is motivated to study hospital 
inpatient flow management and is introduced in [6] (which we will refer to as the “main paper”). To 
analyze this system, a critical step is to obtain the stationary distribution tt for the midnight count 
process {V^., k = 0,1,... }, where denotes the number of customers in the system at the midnight 
of day k, including both customers in service and those waiting in the buffer. To efficiently compute 
TT, especially when the number of servers N is large and the utilization is close to 1, we develop a 
discrete-time diffusion process {X^, k = 0,1,...} to approximate the midnight count process, and 
use the stationary distribution of {X^} to approximate vr. 

In this note, we first prove a limit theorem in Section 1 that supports this diffusion approximation. 
Then, in Section 2 we discuss two methods to numerically compute or approximate the stationary 
distribution of the discrete-time diffusion process {X^}. Finally, in Section 3, we show the accuracy 
of these two methods in approximating tt via numerical experiments. 

1 Diffusion limits for the single-pool model 

Section 4.3 of [6] proposes a discrete-time diffusion process to approximate the midnight count 
process. This approximation is motivated by a limit theorem that shows the convergence of stochastic 
processes. In this section, we prove this limit theorem. 

Instead of fixing the number of servers N, we consider a sequence of Mperi/Geo 2 timeScaie/X 
systems indexed by N, i.e., a sequence of the single-pool models described in the main paper [6]. 
Let be the daily arrival rate of the Xth system. Let m = l//r, the mean LOS, be fixed and 
= {A^m)/N be the traffic intensity of the Xth system. We assume that 

lim A^/X = A*, and lim y/N{l — p^) = /3* for some f3* > 0. (1) 

N^OO N^OO 

Analogous to the conventional many-server queues that model customer call centers [7], we call 
Gondition (1) the Quality- and Efficiency-Driven (QED) condition. 

We use to denote the midnight customer count at the midnight of day k in the Xth system. 
We consider the diffusion-scaled midnight customer count processes = {X^ : k = 0,1,2,...} 


1 


( 2 ) 


for the sequence of singe-pool systems, where for a given k, is defined as 

— 


^/N 

Adapting the derivations in the main paper, we can show that X^ satisfies the following relationship: 


k-l 


yN 






N\- 


)-, A: = 0,l,2..., 


i=0 


( 3 ) 


where 


yN _ -yN I 


^/N 




{pfo,k] ~ + • • ■ + Zk-i)^ + kVN^{p^ — 1 ), 


^(0 A;] ~ -^(0 k] ~ "^^=0 cumulative number of arrivals and departures from 

0 until the midnight (zero hour) of day k in the A^th system, respectively, and = min(XA^,A^) 
is the number of busy servers at the midnight of day i. We assume the initial condition 


X^ ^ X^ as N ^ oo. 


( 4 ) 


where => denotes convergence in distribution. Under the many-server heavy-traffic framework (e.g., 
see [4]), we prove the following limit theorem: 

Theorem 1 Consider a sequence of Mperi/ Geo 2 timeScaie/X single-pool systems that satisfies (1) and 
(4)- For any positive integer K G Z+, X^ ^ X^ on the compact set [0, K] as N ^ oo, i.e., 

(x^,x/^, ...,X^)^ • • • > asN^^. (5) 


The discrete-time limit process X^ = {xj,, A: = 0,1,... } satisfies 

k-l 

4 = + k = 0,l,..., ( 6 ) 

i=0 

where = Y^{k) for A; = 0,1,..., is an embedding of the Brownian motion {Y^{t),t > 0} which 
starts from Xq and has mean —pfi and variance A* + /i(l — p). 

In this limit theorem, we deliberately use the superscript f to differentiate the limit process X^ (and 
the associated process from the diffusion approximation X* (and the associated process Y*) 
introduced in Section 4.3 of the main paper. 

The key step of the proof for Theorem 1 is to show that {Y^,k = 0,1,...} converges to 
{Y^, A: = 0,1,... } on any given compact set [0, K], or equivalently, 

^ (yj,l}t,...,4) asiV^oo. (7) 


Then, the convergence of X^ to X^ naturally follows because of the linear forms in (3) and (6). To 
prove (7), we first prove the convergence of the diffusion-scaled arrival processes in Section 1.1, and 
then the convergence of the discharge processes in Section 1.2. 
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1.1 Arrival process 

For the Nth system, let 
{E^{t), t > 0} defined as 


(o,fc] 


We also introduce a continuous-time process 


E^{t) = ^{E^{t)-K^t), 


( 8 ) 


where E^{■) represents a Poisson process with rate . It is easy to verify that {E^} is an 
embedding of E^{■), i.e., 

E^ = E^{k), fc = 0,l,.... 

Following a standard functional central limit theorem argument, we can show that 

E^{-) => eH-) (9) 

in space D endowed with the Skorohod Ji topology, where E^{-) is a Brownian motion with drift 
0 and variance A*. Because the convergence of stochastic processes implies the convergence of any 
finite-dimensional joint distributions, we then naturally have 


^0 5-^1 5 • • • 5 


eI,eI,...,eI^] as aoo, 


( 10 ) 


where eI = E^{k) is also an embedding of E^{-). 


1.2 Discharge process 

Now we consider the diffusion-scaled discharge processes. For the Ath system, we introduce two 
discrete-time processes: 


D^ = 


1 


y/N 


D 


N 

{0,fc] 


“ + • • • + ) , k — 0,1,. 


N 


and 


= 


\fN 


X] (Ci - a) , = 0,1,..., 

i=l 


where {^j} is a sequence of iid Bernoulli random variables with success probability //. Recall that in 
Appendix C of the main paper, we establish a revised system which tosses coins for every customer 
in service at the midnight to determine the departures each day, and we have proved this revised 
system is equivalent to the original system in distribution. Using the revised system, we can show 
the above two discrete-time processes are equal in distribution, i.e.. 

Thus, it is sufficient to prove for any given K G Z_|_, 


r)N fyN jy 

■^0 ■> 1 ■ ■ ■ 1^ 


{SI,SI,...,S*k) as A^oo. 


( 11 ) 


Here, 5^ = S*{k) is an embedding of the Brownian motion S*{-) with drift 0 and variance ^(1 — /x). 

Let rji = — n, and {rn} forms a sequence of iid random variables with mean 0 and variance 

/u(l — b)- We also define 


fc-i 

i=o 




tA = 




A ’ 


3 



and 


Sn = '^r]i- 


i=l 


Then, we can further rewrite as 


= 


E 


v/iv 

Correspondingly, proving (11) is equivalent to showing 


StNn- 


-^SrpN^, -^SrpN]^, -^SrpNjj^ => (, . . . , S^) US iV — > OO . 
To prove (12), we introduce a continuous-time process {S^{t),t > 0}, where 


( 12 ) 


= t>0. 

In other words, S^{-) is a composition of two continuous processes, and T|j. It is easy to 

verify that {D^} is an embedding of S^{-) because when t = k, T^N = is always an integer 
and 


D^ = 


1 




S* 


If we can show 

in space D endowed with the Skorohod Ji topology as well as 


T^-\ in probability 


(13) 


(14) 


with Tpj = [tj, then applying the random time change theorem, we can prove (12). The convergence 
in (13) follows from the Donsker’s theorem, and we focus on proving (14) below. It is sufficient to 
show for each 0 < k < K, /N —?• 1 almost surely, which we prove with induction. 

We hrst rewrite the system equation under the fluid scaling: 


fc-i 


vN _ 




N\- 


(15) 


i=0 


where 


and 


X{^ -N 


Y'N _ 


N 


m , Kp^ - 1 ) 


+ 


N 


N 


+ 


Vn 


Assume that Xq = N, then Xq = 0 and Zq = N (so Xq —)• 0 is trivial). 
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• When /c = 1, we have 


Recall that {A^ — A^) and rji are centered random variables with mean 0. By the Law of 
Large Numbers, it is obvious that 

—)• 0 a.s. when N ^ oo. 

Thus, —)• 0 a.s., and /N —)• 1 a.s.. 

Assume that at k, we have for all 0 < j < k, XA —)• 0 a.s. and ZAJN —>■ 1 a.s.. Then for 
A; + 1, we have —)■ (fe + 1) a.s. and 


Y 


N 

k-\-l 


i>jv , - A") eS‘ m + i)(p« - 1) 

“ + N + Tv 

E&i(Af 1 - A«) eS‘ m 1 ^ + i)(p« -1) 


N 


rpN 
-‘fc + 1 


N 


y/N 


0 a.s.. 


(16) 

(17) 


Then 

= + a.s., 

j=0 

which completes the proof of Theorem 1. 


2 Computing the stationary distribution of the discrete-time diffu¬ 
sion process 

Motivated by the limit theorem proved in Section 1, Section 4.3 of the main paper [6] proposes a 
discrete-time diffusion process {W^, /c = 0,1,... } to approximate the original midnight count process 
{Xk, /c = 0,1,... }. The dynamics of this approximation process follows: 

fc-i 

Xl=Y,*+fiY.iX:)-, A: = 0,1,2,..., (18) 

i=0 

where Y^ = Y*{k) for A: = 0,1, 2,..., and {Y*(t),t > 0} is a Brownian motion with mean 

9n = A- Nn = p) (19) 

and variance 

a% = A +pNp{l - p) = pNp{2 - p). (20) 

Note that (19) and (20) are different from the mean —pP and variance A* -|- p(l — p) in Theorem 1, 
for two reasons: first, the process X* and Y* are diffusion approximations instead of the limiting 
processes stated in Theorem 1, which is why the term p appears in (19) and (20); second, the 
process X* is to approximate the centered midnight count process (defined as Xk = X^ — N), not 
the diffusion-scaled version as in (2). 

In the next three subsections, we first specify the basic adjoint relationship (BAR) for this 
discrete-time diffusion process X*. Then, we discuss two ways to numerically calculate/approximate 
the stationary distribution of X*: (i) a projection algorithm that numerically solves the BAR, and 
(ii) an approximate formula. 
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2.1 Basic adjoint relationship 

The state space of fe = 0,1,... } is M. One can check that {X'^, /c = 0,1,... } is a discrete-time 
Markov process, since 


- XI = nVi - + ^l{Xl)-, for fe = 0,1,..., 


and ~^k : ^ = 0,1,...} is a sequence of iid normal r.v. with mean 07 v and variance The 

transition density of the Markov process is 


p{x,y)=¥{Xl^, = y\Xl = x) 


when X > 0, 
(y - (1 - y)x ), when x < 0, 


( 21 ) 


where (t>g denotes the normal density function with mean 9 and variance Let C'b(K) denote 
the set of bounded, continuous functions on M. For each / G C;,(M), define 


P/(x) 


p{x,y)f{y)dy 


for each x G M. 


One can check that P/ G C'b(M). It follows that the stationary density 7r(x) satisfies 


P/(x)7r(x)(ix = / f{x)TT{x)dx, V/ G 


( 22 ) 


or equivalently. 


L/(x)7r(x)(ix = 0, 


V/ G Cb{R), 


(23) 


with L/(x) = P/(x) — /(x). We call (23) the basic adjoint relationship (BAR) that governs the 
stationary density of the discrete-time Markov process {W^, A: = 0,1,... }. 


2.2 A projection algorithm 

The BAR (23) is in the same format as (2.5) of [5]; the latter BAR is for the stationary density 
of a (continuous-time) diffusion process. As such the algorithm developed in [5] can be applied to 
compute the stationary density vr* of the discrete-time diffusion process {X^,k = 0,1,...}. We 
outline the algorithm here, commenting on the differences when appropriate. 


2.2.1 Reference density and the space L^(M, r) 

To compute the stationary density tt* , we first need a reference density r such that 

/ r{x)dx = 1. 

Jr 

We use the approximate formula tt in Section 4.3.2 of the main paper (also see 34 below) as the 
reference density r. 

Next, we define the ratio function as: 

q{x) = ^ for X G M. (24) 

r[x) 

With the given reference density r, if we can compute the ratio function q, then we can compute 
the stationary density via 

7r*(x) = q{x)r{x) for x G M. 
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To compute g, we plug (24) into (23) and get 


Ijf {x)q{x)r{x)dx = 0, V/ G C'fe(M). 


(25) 


Following the notation in [5], we use L^(M, r) to denote the space of all square-integrable functions 
on M with respect to the measure that has density r. Namely, L^(M, r) is the set of measurable 
functions / on M that satisfy 

/ f‘^{x)r{x)dx < oo. 


We adopt the same inner product on L^(]R, r) as in [5], that is, 

(/,/) = / f{x)f{x)r{x)dx, for fj G L‘^{R,r). 
Jr 


(26) 


In (3.2) of [5], the authors made an important assumption on the reference density. Namely, they 
assumed that the reference density was chosen so that 


q G L^(M, r). 


(27) 


With our choice of the reference density r, we have been unable to verify that condition (27) is 
satisfied. We leave it as a conjecture that condition (27) is satisfied. The remainder of this section 
assumes that the conjecture is true. 


2.2.2 Orthogonal projection 

Note that the BAR (25) is equivalent to 

(L/, q) = 0 for each / G C;,(]R). 


Thus, q satisfying the BAR is equivalent to q being orthogonal to L/ for each / G C'b(M). We define 
a space H as 

H = the closure of {L/ : / G C'ft(M)}, 

which is a subspace of L^(M, r). Therefore, q satisfying the BAR is equivalent to q being orthogonal 
to space H. Therefore, our task is to find a function q that is orthogonal to space H. To do so, we 
consider a constant function e with e{x) = 1 for each x G M. Since 

{e,q) = / e{x)q{x)r{x)dx = / 'K*{x)dx = l, (28) 

1/ M. M. 


one can check that e ^ H because otherwise (e,g) = 0, contradicting (28). We use e to denote the 
projection of e onto H. Then, e — e 7 ^ 0 and it must be orthogonal to H. Once we have e, we obtain 
the ratio function q by 

e — e 


q = 


e — 


2 ’ 


where || • || is the induced norm from the inner product (26) with ||/|p = (/, /) for / G L^(M, r). 
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2.2.3 Finite-dimensional approximation 

The projection of e onto H can be expressed as 

e = argmin^jgji^lle -/i||. (29) 

The space H is linear and inhnitely dimensional. To compute the projection numerically, we use a 
hnite-dimensional subspace to approximate H and hnd the projection of e on namely, 

Cfc = argmin^jgj:^^I|e - ^11• (30) 

Let Ck be a hnite-dimensional, linear subspace of C;,(]R). Then Hj. = {L/ : / G is a hnite- 
dimensional subspace of H. Assume that {fi : i = 1,2,... ,m} G Ck is a basis of Ck- Then, since 
the projection Bk G Hk, it can be represented as a linear combination of {L/j : i = 1,2,... ,m}. 
That is, 

m 

ek = '^aiLfi (31) 

i=l 

where G M for i = 1,2,... ,m. 

To compute the vector of coefficients a = (ai,..., am)^ we use the fact that (e — e^, L/j) = 0 
for i = 1, 2,..., m. Consequently, we obtain a system of linear equations 

Aa = 13, (32) 

where Aij = (L/j, L/j) and /3j = (e, L/j) for i,j = 1,..., m. The matrix is symmetric, semi-positive 
dehnite, but can be singular. Although the solution to the system of linear equations may not be 
unique, projection Bk is unique. When A is singular or nearly singular, one can solve (32) by direct 
methods such as the QR decomposition and the Cholesky decomposition or by iterative methods 
such as LSQR [10]. The Cholesky decomposition exploits the symmetric and semi-positive dehnite 
properties of A even when A is singular [1, 8], whereas QR decomposition does not. Unlike many 
other iterative methods, LSQR can handle matrix A when it is singular. LSQR does not exploit 
semi-positive dehniteness. 

Once we get the vector of coefficients a = (ai,..., am)' by solving the system of linear equations 
(32), we can compute Bk as in (31). Eventually, we can approximately compute the stationary 
density tt* as 

7r*(x) r(x)(i^—Vx G M. (33) 

e - er b 


2.2.4 FEM implementation 

In our implementation, we use the hnite element method (FEM) to construct the approximate space 
Ck, following Section 3.3 of [5]. The numerical results in this paper for approximating the stationary 
density vr with the projection algorithm all follow this FEM implementation. In Proposition 3 of 
Dai and He [5], they proved the convergence of using (33) to approximate vr as Hk t H. Their proof 
applies to our setting when (27) is satished. 


2.3 Approximate formula for the stationary density 


In Section 4.3 of the main paper, the following formula vf is proposed as a proxy for the stationary 
density vr* of the diffusion process X*\ 


vf (x) 


ai exp(20jva;/(T^), x > 0; 

a 2 exp(-( 2 // - iJ?){x - 6 n// 2a%), x < 0; 


(34) 
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where ai and 02 are normalizing constants that make ^{x) continuous at zero and j’jg7r(a:)dx = 1. 

As mentioned in the main paper, the rationale of this approximate formula is based on the 
analogy between : /c = 0,1, 2 ...} and > 0}, where 

X{t) = Y{t) +n [ {X{s))-ds, t>0, (35) 

Jo 

and {Y{t),t > 0} is a Brownian motion. To get the stationary density of X, Browne and Whitt [3] 
have suggested that since (i) A is a Ornstein-Uhlenbeck (OU) process on (— 00 ,0] and the stationary 
density of an OU process has a Gaussian form and (ii) A is a reflected Brownian motion (RBM) on 
[0, 00 ) and the stationary density of a RBM has an exponential form, then the stationary density of 
A can be obtained by piecing together the Gaussian and exponential densities. 

We use the same piecing technique in our setting. Specihcally, {A^ : k = 0,1,2...,} behaves 
as a discrete version of the OU process on (— 00 , 0 ] and as a reflected random walk on [ 0 , 00 ). We 
show in Proposition 1 below that the stationary density of the discrete-time OU (DOU) process 
also has a Gaussian form. For the reflected random walk, existing research shows that it has an 
exponential tail [9, 11, 2]. Therefore, we piece together a Gaussian density and an exponential 
density and propose using (34) to approximate vr*. In the next two subsections, we hrst prove that 
the stationary density of the discrete-time OU process has a Gaussian form, then we show the details 
of deriving formula (34). 


2.3.1 The stationary distribution of a discrete OU process 

Similar to the continuous-time version of the Ornstein-Uhlenbeck process, we dehne its discrete-time 
version {ADOU^fc = 0,l,...} as: 

ADOU=yDOU_^^^DOU^ A: = 0,1,... (36) 

i=0 

where := ^ = 0) !> ■ ■ • } is a Gaussian random walk, i.e., {^i} is a sequence of iid 

random variables following a normal distribution with mean 9 and variance 

The following proposition says the stationary density for a discrete OU process has the Gaussian 
form, which is consistent with that in a continuous-time OU process. 

Proposition 1 Given 0 < /x < 1, for a discrete-time Ornstein-Uhlenbeck (DOU) process {A^®^, k = 
0,1,... } satisfying 

X^OU A: = 0,1,... (37) 

i=0 

where {YjP^^^} is a Gaussian random walk with drift 6 and variance the stationary density of 

2 

the DOU process, tt, is a normal density with mean 9/fi and variance • 

Proof for Proposition 1. Note that the DOU process {A^®^,A; = 0,1,...} satisfying (37) is a 
Markov process since 

V-DOU vDOU _ taDOUn ,, vDOU 

~ — 1-'fc+l “ -'fc ) ~ 

The transition probability from state y to state x is 

1P(A^°X = x\X^^^ = y) = - (1 - iM)y), 
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where (j)g (s) denotes the probability density function associated with a normal random variable 
with mean 9 and variance 

To prove this proposition, we just need to show that 


7r(x) 



F{x\y)7r{y)dy 


(38) 


for any given x, where 


y?) 

TT[x) = —^==-• exp 


\/^i 


vrcr 


{2y - y^){x - 9/yf 


We have 

'^{x\y)TT{y) = 


^J2y-y?' 1 

\/27r cr \f^<J 
^[2y-y^) 1 

^[2y-y^) 1 


exp 

exp 

exp I - 


{2y - y'^){y - e!yf 
2a2 


2ct2 


exp 


y'^ - 2[{l - y)x + e]y 
2a2 

[y-{{l-^,)x + e)Y 

2a2 


exp 


{x- (I- y)y-ef 
2a2 

{2y — y^)9'^/y^ + (x — 0)^ 
2(^2 


• exp - 


(2y - y^)9^/y^ + (x - 6y - [(1 - y)x + 0]" 


Among which, 


V (x) = exp ( — 


{2y - y^)9‘^/y^ + (x - 9f - [(1 - y)x + 0]^ 


= exp 


= exp 


2ct2 

{2y — y^)9‘^ / y^ + {2y — y?)x‘^ — 2(2 — y)9x 


2ct2 


{2y - y^){x - 9/yf 
2ct2 


Then, we have 


[x\y)'K[y)dy = ^—=-exp 

y l7:(7 

■f 

J — oo V 

\/(2/i-//2) 




TTCT 


• exp 


{2y - y'^){x - 9/yf 
2ct2 

2^2 

{2y-y‘^){x-9lyf 

2c72 




which takes the exact form as the normal density with mean 9/y and variance l(2y — y^) and 
thus, equals to 7r(x). This completes our proof for vr being the stationary density. 


2.3.2 Derivation of the approximate formula 

Based on Proposition 1, we conjecture that the stationary distribution of X* can be approximated 
by the following form: 


if (x) 


7ri(x) = ai exp(—7x), x > 0; 

7r2(x) = a 2 exp(-^?^hh^^^±^)^ x < 0. 


(39) 
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For the ease of exposition, we use 6 and cj^ instead of 6 n and to denote the mean and variance of 
the discrete-time diffusion process X*. Moreover, in (39), ai and 02 are two normalizing constants, 
7 is the unknown parameter for the exponential density part, and we define 

/3 = 


The stationary density should satisfy 

/ OO 

p{y, x)7r{y)dy, 

-OO 

one special form of the BAR (22), or equivalently, 

poo pO 

Tf{x)= p{y,x)TTi{y)dy + p{y, x)Ti 2 {y)dy, 

7o J—OO 


(40) 


where p{y,x) is the transitional density of X* (from state y to state x) defined in (21). 
We rewrite Equation(40) as follows. First, for y > 0, we have 


p{y,x)TTi{y) = 


ai 


\/^a 

ai 

cti 

cti 

OL\ 


\/^i 


• exp 

• exp 

• exp 

• exp 

• exp 


{x-y^ 

2ct2 


• exp(-7y) 


- 2(x + - cr'^'y)y (x -h y/df 


na 


2ct2 

[y-{x + fj,/d- 
2a^ 

[y-{x + p/3- cj 27)]2 
2c72 

. . exp 


• exp 

• exp 


(x -k pfdy - (x + p/d - 0 -^ 7 )^ 
2^2 

cr27[2x -I- {2p(d - 0 - 27 )] 


2ct2 

[y-{x + pfd- 
2ct2 


Therefore, 


poo 

/ p{y,x)'Ki{y)dy = ai 
Jo 


exp 


= ai exp 
Second, for y < 0, we have 


7(2y/3 - 0-27) 


7(2y/3 - 0-27) 


exp (— 7 x). 

[y-{x + p/d- 

2cj2 


^ exp(- 7 x) • [1 - (-X - (2/i/3 - 0 -^ 7 ))] 


p{y,x)'K2{y) = 


Therefore, 


a2 


\[^o 

0.2 


exp 




'KO 


■ exp 


(x - (1 - y)y-ky/3)^ 
2a2‘ ^ 

(y - ((1 - y)x -/i/3))2 
2ct2 


exp 


• exp 


(2y-y2)(y + ;5)^ 


2ct2 

(2/x-y2)(^ + ;3)2 

2c72 


rO 


p{y,x)'K 2 {y)dy = 


OL2 


I—OO ^TTi 


exp 


Tra 


(y - ((1 - y)x -/i/3))^ 


2(72 


exp 


= 02 exp 

= 02 exp 


{2p-p?){x^fdf 

2(t2 

(2y-y2)(^ + ;5)2 

2(t2 


/■u ^ 

/ 

J—OO V "Trcj 
) • ^-t,fi,ak-X + PX). 


(2y-y2)(x + /3)2 
2(t2 

(y- {1- p)x + pld)^ 
2(72 


dy 

dy 
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If (40) holds, when x > 0, we should have 

7(2^/3 - 0 - 27 ) 


aiexp(- 7 x) = ai exp 7) ^ exp(- 7 x) • [l - (-x - (2^/3 - q-%))] 

/ (2^ —/i2)(x +/3)2\ 

+ 02 exp (-- •$_^^_^ 2 (-x + //x), 


which is equivalent to 

«!exp(— 7 x) 

(2m-m')(x + /3)2 


2^x2 


1 - exp ( I (1 _ - (2/i/3 - 0-%))) 


= a 2 exp 


2a2 




(41) 


Similarly, if (40) holds, when x < 0, we should have 


a 2 exp 


(2/i-^2)(3,_^/5)2 

2^2 


= a\ exp 


7(2^/3 - 0 - 27 ) 


^ exp(- 7 x) • [1 - 4 >-;,/ 3 ,a 2 (-x - (2ir/3 - ct 


(2/r —/i2)(x +/3)2' 

+ 02 exp ( -(-X + /xx). 


2X72 


which is equivalent to 
ai exp 

= a 2 exp 


7(2/x/3 - (727) 


{2ii - ij?){x + (5) 

2X72 


^ exp(- 7 x) • [1 - 4'_^^_^2(-x - (2/X/3 - X 727 ))] 

• (1 - + nx)) . 


«!• 


When X = 0, Equations (41) and (42) become 
7(2/x/3 - X727)' 


1 — exp — 


(1 - - a^j))) 


= 02 exp 


(2^-/x2)/32 
2X72 


(42) 


(43) 


and 
a I exp 


7(2/x/3 - X727) 


y [1 - 4>_^,g,^2(-(2^/3-X727))] = 02 exp 2 x 7 ^ ^^ )•( 


1 _ $ 


— fj,l3,a-^ (0)) I 

(44) 


respectively. 

Recall that 7r(x) is continuous at x = 0. Thus, the two normalizing constants satisfy: 

7ri(0) = 01 = 02 exp ^ = 7 x 2 ( 0 ). (45) 

Comparing (45) with (43) and (44), we find that when 

X727 = 2 n (3 = - 26 , 

or equivalently, 

26 , , 

7 = -^, (46) 

both (43) and (44) can be satisfied. Therefore, we choose 7 in (46), which eventually gives us (34). 
Unfortunately, using this 7 , we are unable to show (41) and (42) hold for a general x. 


^))] 
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(a) N = 500, A = 90.95, p = 0.96 



(b) iV = 66, A = 11.37, p = 0.91 



0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 

Midnighi Count 

(c) N = 18, A = 3.03, p = 0.89 


Figure 1: Stationary distribution of the midnight customer count from exact Markov chain analysis and 
diffusion approximations. Here, the mean LOS is 5.3 days, and we do not specify the discharge distribution 
because it does not affect the midnight customer count distribution. 




(a) Eoo(0(t)), N = 500 


(b) Eoo(fF(t)), AT = 500 


(c) Poo{W{t) > 6/24), N = 500 


Figure 2: Time-dependent performance curves from exact analysis and diffusion approximations. Here, 
A = 90.95 for N = 500. We fix the mean LOS as 5.3 days and use the baseline discharge distribution. The 
three performance curves in each subfigure are from normal approximations using (i) tt solved from exact 
Markov chain analysis, (ii) tt in (34), and (iii) tt* solved from the projection algorithm, respectively. 


3 Numerical results on diffusion approximations 

3.1 Approximation for the midnight count distribution 

Figure 1 compares the stationary distributions of the midnight customer count solved (i) from the 
exact Markov chain analysis, (ii) from using the approximate formula fr in (34), and (iii) from 
using the projection algorithm specified in Section 2.2. The parameter settings for these numerical 
experiments are the same as those in Section 5 of the main paper. We test a large system (A = 500) 
and two small systems (A = 66 and 18), with the utilization p being 96%, 91% and 89%, respectively. 

3.2 Time-dependent performance 

Figures 2, 3 and 4 show the time-dependent performance for systems with A = 500, 66, and 18, 
respectively. The three curves in each subfigure are obtained from normal approximations using (i) 
TT solved from exact Markov chain analysis, (ii) fr in (34), and (iii) vr* solved from the projection 
algorithm specified in Section 2.2. 
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(a) Eoo(Q(i)), AT = 66 (b) E^{Wit)), N = 66 (c) P^{W(t) > 6/24), Af = 66 

Figure 3: Time-dependent performance curves from exact analysis and diffusion approximations. Here, 
A = 11.37 for N = 66. We fix the mean LOS as 5.3 days and use the baseline discharge distribution. The 
three performance curves in each subfigure are from normal approximations using (i) tt solved from exact 
Markov chain analysis, (ii) tt in (34), and (iii) tt* solved from the projection algorithm, respectively. 



(a) Eoc(Q(t)), AT = 18 (b) Eoc(W(t)), A^ = 18 (c) P^{W{t) > 6/24), M = 18 


Figure 4: Time-dependent performance curves from exact analysis and diffusion approximations. Here, 
A = 3.03 for N = 18. We fix the mean LOS as 5.3 days and use the baseline discharge distribution. The 
three performance curves in each subfigure are from normal approximations using (i) tt solved from exact 
Markov chain analysis, (ii) tt in (34), and (iii) tt* solved from the projection algorithm, respectively. 
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